{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Binary Classification with Support Vector Machines"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "by Soeren Sonnenburg"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook illustrates how to train a binary support vector machine (SVM) classifier with shogun. \n",
      "A classifier attempts to distinguish objects of different type. In case of of binary classification there are just two types of objects that we want to distinguish."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "More formally, we want to learn a decision function $f(x) \\mapsto \\{+1,-1\\}$ based on a number of training examples $(x,y)_{i=0}^{N-1}$ where $x\\in\\mathcal{X}$ and $y\\in\\{-1,+1\\}$. \n",
      "Now a [SVM](https://en.wikipedia.org/wiki/Support_vector_machine) is a binary classifier that tries to separate objects of different classes by finding a (hyper-)plane such that the margin between the two classes is maximized:\n",
      "\n",
      "More formally, a support vector machine is defined as\n",
      " $$f({\\bf x})=sign\\left(\\sum_{i=0}^{N-1} \\alpha_i k({\\bf x}, {\\bf x_i})+b\\right)$$\n",
      " \n",
      " where $N$ is the number of training examples\n",
      " $\\alpha_i$ are the weights assigned to each training example\n",
      " $k(x,x')$ is some kernel function\n",
      " and $b$ the bias.\n",
      " \n",
      " Using an a-priori choosen kernel, the $\\alpha_i$ and bias are determined\n",
      " by solving the following quadratic program\n",
      " \n",
      "  \\begin{eqnarray*}\n",
      "       \\max_{\\bf \\alpha} && \\sum_{i=0}^{N-1} \\alpha_i - \\sum_{i=0}^{N-1}\\sum_{j=0}^{N-1} \\alpha_i y_i \\alpha_j y_j  k({\\bf x_i}, {\\bf x_j})\\\\\n",
      "       \\mbox{s.t.} && 0\\leq\\alpha_i\\leq C\\\\\n",
      "                   && \\sum_{i=0}^{N-1} \\alpha_i y_i=0\\\\\n",
      "\\end{eqnarray*}\n",
      " where C is a pre-specified regularization parameter trading of the size of the margin and the training error.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let us see how one can train a support vector machine with shogun:\n",
      "\n",
      "In a first step we just generate some two-dimensional Gaussians.\n",
      "\n",
      "$x_-\\sim{\\cal N_2}(0,1)-d$\n",
      "\n",
      "$x_+\\sim{\\cal N_2}(0,1)+d$\n",
      "\n",
      "and corresponding positive and negative labels. We create traindata and testdata with ```num``` of them being negatively and positively labelled in traindata,trainlab and testdata, testlab. For that we utilize shoguns gaussian mixture model class (GMM) from which we sample the data points and plot them."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from modshogun import *\n",
      "\n",
      "num=1000;\n",
      "dist=1.0;\n",
      "\n",
      "gmm=GMM(2)\n",
      "gmm.set_nth_mean(array([-dist,-dist]),0)\n",
      "gmm.set_nth_mean(array([dist,dist]),1)\n",
      "gmm.set_nth_cov(array([[1.0,0.0],[0.0,1.0]]),0)\n",
      "gmm.set_nth_cov(array([[1.0,0.0],[0.0,1.0]]),1)\n",
      "\n",
      "gmm.set_coef(array([1.0,0.0]))\n",
      "xntr=array([gmm.sample() for i in xrange(num)]).T\n",
      "xnte=array([gmm.sample() for i in xrange(num)]).T\n",
      "\n",
      "gmm.set_coef(array([0.0,1.0]))\n",
      "xptr=array([gmm.sample() for i in xrange(num)]).T\n",
      "xpte=array([gmm.sample() for i in xrange(num)]).T\n",
      "\n",
      "\n",
      "traindata=concatenate((xntr,xptr), axis=1)\n",
      "testdata=concatenate((xnte,xpte), axis=1)\n",
      "\n",
      "trainlab=concatenate((-ones(num), ones(num)))\n",
      "testlab=concatenate((-ones(num), ones(num)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "_=scatter(traindata[0,:], traindata[1,:], c=trainlab, s=100)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using this data we create shogun features and label objects:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "feats_train=RealFeatures(traindata)\n",
      "labels=BinaryLabels(trainlab)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Just for fun we compute the kernel matrix and display it - the nice block structure of the matrix suggests that the data will be nicely separable (examples 0..999 when compared to each other have mostly higher scores than when compared to example 1000...1999 and vice versa)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "width=2\n",
      "kernel=GaussianKernel(feats_train, feats_train, width)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "km=kernel.get_kernel_matrix()\n",
      "_=imshow(km)\n",
      "_=colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we train an SVM with this GaussianKernel. We use LibSVM but we could use any of the [other SVM](http://www.shogun-toolbox.org/doc/en/current/classshogun_1_1CSVM.html) from shogun. They all utilize the same kernel framework and so are drop-in replacements."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C=1.0\n",
      "svm=LibSVM(C, kernel, labels)\n",
      "_=svm.train()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We could now check a number of properties like how many support vectors the trained SVM has, i.e. how many of the $\\alpha_i$ above are non-zero:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print svm.get_num_support_vectors()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "or what the value of the objective function returned by the particular SVM learning algorithm or the explictly computed primal and dual objective function is"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "libsvm_obj=svm.get_objective()\n",
      "primal_obj, dual_obj=svm.compute_svm_primal_objective(), svm.compute_svm_dual_objective()\n",
      "\n",
      "print libsvm_obj, primal_obj, dual_obj"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "and based on the objectives we can compute the duality gap, a measure of convergence quality of the svm training algorithm. In theory it is 0 at the optimum and in reality at least close to 0. The"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"duality_gap\", dual_obj-primal_obj"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's now compute the test error"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "out=svm.apply(RealFeatures(testdata))\n",
      "\n",
      "evaluator=ErrorRateMeasure()\n",
      "print \"Test error is %2.2f%%\" % (100*evaluator.evaluate(out,BinaryLabels(testlab)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now lets plot the contour output on a $-5...+5$ grid for \n",
      "\n",
      "1. The Support Vector Machines decision function $\\mbox{sign}(f(x))$\n",
      "2. The Support Vector Machines raw output $f(x)$\n",
      "3. The Original Gaussian Mixture Model Distribution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size=100\n",
      "x1=linspace(-5, 5, size)\n",
      "x2=linspace(-5, 5, size)\n",
      "x, y=meshgrid(x1, x2)\n",
      "grid=RealFeatures(array((ravel(x), ravel(y))))\n",
      "grid_out=svm.apply(grid)\n",
      "z=grid_out.get_labels().reshape((size, size))\n",
      "\n",
      "figure(figsize=(20,5))\n",
      "subplot(131)\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n",
      "z=grid_out.get_values().reshape((size, size))\n",
      "\n",
      "subplot(132)\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)\n",
      "\n",
      "subplot(133)\n",
      "gmm.set_coef(array([1.0,0.0]))\n",
      "gmm.set_features(grid)\n",
      "grid_out=gmm.get_likelihood_for_all_examples()\n",
      "zn=grid_out.reshape((size, size))\n",
      "gmm.set_coef(array([0.0,1.0]))\n",
      "grid_out=gmm.get_likelihood_for_all_examples()\n",
      "zp=grid_out.reshape((size, size))\n",
      "z=zp-zn\n",
      "c=pcolor(x, y, z)\n",
      "_=contour(x, y, z, linewidths=1, colors='black', hold=True)\n",
      "_=colorbar(c)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And voila! The SVM decision rule reasonably distinguishes the red from the blue points. Despite being optimized for learning the discriminative function maximizing the margin, the SVM output quality wise remotely resembles the original distribution of the gaussian mixture model."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}